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♦ 1. Introduction: 



We present a software package DiracQ , for use in quantum many-body 
Physics. It is designed for helping with typical algebraic manipulations that 
arise in quantum Condensed Matter Physics and Nuclear Physics problems, 
and also in some subareas of Chemistry. DiracQ is invoked within a Math- 
ematica session, and extends the symbolic capabilities of Mathematica by 
building in standard commutation and anticommutation rules for several 
objects relevant in many-body Physics. It enables the user to carry out com- 
putations such as evaluating the commutators of arbitrary combinations of 
spin, Bose and Fermi operators defined on a discrete lattice, or the position 
and momentum operators in the continuum. Some examples from popular 
systems, such as the Hubbard model, are provided to illustrate the capabili- 
ties of the package. 

A word about the underlying philosophy of this program is appropriate 
here. Physicists approach calculations in quantum theory with a certain in- 
formality that contrasts with the approach of mathematicians, who usually 
insist on a somewhat more rigorous (notational) framework. Some formality 
is also present in the method of programming a computation with older com- 
puter languages such as Fortran. The physicist's informality refers specifi- 
cally to a deferment of definitions and properties of objects, until one actually 
needs them. This enables one to write down and manipulate and to com- 
pound expressions, often into simple and useful final results. In our view the 
physicist's approach owes much to the notation introduced by Dirac in 1927. 
In the classic paper "The Physical Interpretation of Quantum Dynamics", 
P. A. M. Dirac Proc. Roy. Soc.A113,621(1927), we find the first mention of 
c-numbers and q-numbers as follows: 

. . . one cannot suppose the dynamical variables to be ordinary numbers 
(c-numbers), but may call them numbers of a special type (q-numbers). 

This inspired notation helps us to treat the commuting parts of the ex- 
pression with standard care (due their c-number nature), while reserving 
extra care for handling the q-number parts. 

In DiracQ , we initially declare a set of symbols to be operators, i.e. q- 
numbers. The properties of standard operators such as Bose, Fermi, position- 
momentum, angular momentum etc are already programmed, and if required, 
a user could define more exotic operators with specific properties. Once this 
is done, any expression is split into its c-number and q-number parts as 
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the first and fundamental operation. Complicated input expressions may 
be written with a fairly informal syntax, with mixtures of c-numbers and q- 
numbers, and as the sums of such expressions. DiracQ handles them by first 
separating them into their c-number and q-number parts. Standard opera- 
tions on operators, such as adding, multiplying or commuting them is then 
straightforward, the c-number parts are spectators for most part while the 
q-numbers are combined using the given rules. Finally the expressions are 
written back in standard input like notation. The package DiracQ consists 
of a set of mutually dependent functions to perform most of the operations. 
These functions are most often named using the last letter as "Q" ; for exam- 
ple we denote a function (described later) as SimplifyQ[a], thus distinguishing 
it from the function Simplifyfa] of Mathematica . 

This version of DiracQ uses commands that call upon the symbolic capa- 
bilities of Mathematica, and exploits its ability to deal with non commuting 
symbols. It might be possible to couple the commands of DiracQ with an 
OSS program such as Sage, in view of the relatively small fraction of com- 
mands of Mathematica that are actually used here. The present version 
of DiracQ was developed and tested on Version 8 of Mathematica, and re- 
quires the palettes feature of this version in order to declare the standard 
non commuting symbols. A version avoiding the palettes is planned for the 
future. 

♦ 2. Contents: 

The package distribution contains the following files: (X is the version 
number) 

1. Introduction.pdf 

This article. 

2. DiracQ.VX.m 

The main package implementing DiracQ . 

3. GettingStarted VX.nb 

A first introduction to loading DiracQ and simple examples. 

4. Glossary _VX.nb 
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A list of all commands in DiracQ , their description and a simple 
example or two of their usage. It also contains, under the heading 
DiracQPalette, a description of the palette used in this notebook, 
with its various predefined operators and instructions on how to "turn 
them on" . 

5. Tutorial VX.nb 

A more detailed introduction to the commands and their illustration in 
simple problems. This notebook has some overlap with the notebook 
Glossary _VX.nb, it provides an in depth commentary on the usage of 
the various functions defined in the package. Simple examples with 
harmonic oscillator operators are provided as warming up. Examples 
are provided for all the functions available in DiracQ . Users who 
would like to enlarge the class of operators and define their own set, 
would find a helpful example treated, where the defining relations of 
the Virasoro algebra are added, and a few simple computations with 
their generators carried out. 

6. Examplebook_VX.nb 

A set of calculations using of DiracQ in problems of various level of 
difficulty Harmonic oscillator commutations provide a simple introduc- 
tion to the package, followed by a demonstration of the conservation 
of the Runge-Lenz-Laplace vector in the Hydrogen atom. The Bra and 
Ket objects of DiracQ are introduced. Fermionic operators are illus- 
trated in the context of the Hubbard model, where the commutators of 
combinations of Fermi operators with the Hamiltonian are evaluated, 
leading to exact expressions for the exact second and third moments of 
the electron spectral function. Applications to spin half particles and 
Pauli matrices are given using R. J. Baxter's celebrated proof (1971) 
of the integrability of the 8-vertex model. The relations between the 
Boltzmann weights required for satisfying the star triangle relations 
are recovered using DiracQ . Also provided is a more complex exam- 
ple based on the 1-d Hubbard model. Its S matrix ( Shastry 1986) 
is known, which can further be used to construct an inhomogeneous 
generalization of the Hubbard model. This construction is dependent 
upon the S matrix satisfying a more complex relation that is very hard 
to check analytically (due to the large number of terms involved). It 



4 



can be verified easily using DiracQ , providing a nontrivial demonstra- 
tion. Another example illustrates the higher conservation laws of the 
1-d Hubbard model that are non quadratic in the Fermi operators. 

DiracQ can also be used in the initial stages of numerical calculations 
on finite systems, since it can generate the numerical Hamiltonian ma- 
trix on a given cluster. This is illustrated in a 4 site cluster with nearest 
and next nearest hopping bonds in the sector of 2 up and 2 down parti- 
cles. This example also illustrates the strategies for generation of states 
using the Bra and Ket objects of DiracQ . 



♦ 3. Contacts: 

The authors can be contacted with comments and suggestions as well as 
further user generated examples of DiracQ by many channels. The preferred 
one is the email address for this purpose : 

DiracQ@gmail.com 



♦ 4. Suggested Acknowledgement of DiracQ 

If DiracQ is useful in obtaining results leading to publication, a citation 
would be appropriate and appreciated. 

Suggested Citation : 

This work contains results obtained by using DiracQ: A Quantum Many- 
Body Physics Package , J. G. Wright and B. S. Shastry, arXiv:1301. (2013). 



♦ 5. Acknowledgement: 

This project was supported in its preliminary stage by DOE under Grant 
No. FG02-06ER46319. 



5 



DiracQ Example Book 



A Compilation of Problems in Quantum Theory Performed Using the DiracQ Package 

The following problems demonstrate how the DiracQ package can be used to solve non-trivial problems. Throughout 
this notebook ensure that you have the relevant operators activated on the DiracQ palette before evaluating any input. 

SetDirectory [NotebookDirectory [ ] ] ; 
Get [ "DiracQ_VO.m" ] ; 

■ (I) Canonical Pairs 

The package contains canonical position and momentum operators. Their use is demonstrated below for some simple 
examples. 

?P 



p is the canonical momentum operator. This operator can be called with one argument, taken to be 
site index, or two arguments. The second argument will be taken to be coordinate direction. Also 
included is the 3 dimensional canonical momentum vector, represented by OverVector[p], or p. 

?q 



q is the canonical position operator. This operator can be called with one argument, taken to be 
site index, or two arguments. The second argument will be taken to be coordinate direction. 
Also included is the 3 dimensional canonical position vector, represented by OverVector[q], or q. 

The physical system can be in arbitrary dimensions and may refer to an arbotrary number of particles. The canonical 
pairs are endowed with indices: either one or two indices can be used. The first index is taken to be site label. The 
second index is taken to be the component. If only one index is given, a 1 -dimensional system is assumed. This is 
demonstrated below. [Note: unless you choose to declare q and p as active operators in the DiracQ Palette by ticking 
the boxes in the Palette, (or more aggressively turn all symbols on), the following commutators will turn out to vanish. 
If that happens, go back and turn them "on". ] 
Commutator [p [j ] , q[i]] 

- i ft 6 [if j ] 

Commutator [p[ j , a], q[i, 13] ] 

-ifi6[i, j] <5[a, /3] 
Commutator [p[ j , y] , q[i, y] ] 
-ih 5[i, j ] 
Commutator [p, q] 



Comment: In the last output, notice that p and q without any indices are treated as c-numbers. This can be a slight 
nuisance, but one resolution is to treat a set of harmonic oscillators indexed by a "i". 
Clear [i, j ] 

p[i] 2 mw 2 

H = + q[i] 2 ; 

2 m 2 



2 | Examplebook_VO.nb 



Commutator [H, p[i]] 

i (mo) 2 fij ** q [i] 

Next we construct creation and destruction operators from the canonical pairs. 

P[i] \ 



a[i_] 



m 01 
2 h 



m oi 

at[i_] := A / 

2 h 



q[i] + i 



q[i] - i ■ 



m oi ) 

P[i] 
m u 



Commutator [a [i] , at[i]] 

1 

n[i_] : = at[i] **a[i] 

Here we show that [n[i], a[i]] = -a[i] and similarly [n[i],at[i]]=at[i]. SimplifyQ puts the expressions in a standard 
form where the identity of two expressions becomes clear. 
SimplifyQ [a [i] ] 




P[i] 



** qLiJ 



'2 V2 
SimplifyQ [Commutator [n [i] , a[i]] +a[i]] 



SimplifyQ [Commutator [n [i] , at [i] ] - at [i] ] 



(II) Hydrogen atom and Runge Lenz Laplace vectors 

Here we define the Hamiltonian of the Hydrogen atom. We demonstrate the conservation of the Runge-Lenz-Laplace 
vector by using the package. 

We use the three dimensional vector notation for the canonical pair. For example 
q[i] 

{q[i, x] , q[i, y] , q[i, z] } 

q[i] -q[i] 

q[i, x] 2 +q[i, y] 2 +q[i, z] 2 
Clear [H] ; 



p[i].p[i] 

H = Z 

2 m 

? NCcross 



-q[i] 



Non Commutative cross product of two 3dimensional vectors retaining the order of the operators. 

L[i_] = NCcross [q [i] , p[i]] 

{q[i, y] **p[i, z] -q[i, z] **p[i, y] , 
-q[i, x] **p[i, z] +q[i, z] **p[i, x] , q[i, x] **p[i, y] -q[i, y] **p[i, x] } 
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Commutator [£,[i] [ [1] ] , h] 







Commutator [£[i] [ [2] ] , h] 







SimplifyQ[commutator[L[i] [ [3] ] , h] ] 







The answer is so simple that it looks as though we have not performed any operation whatsoever. However, if we 
leave basic commutators unevaluated we see that indeed the computation was not entirely trivial (feel free to try it by 
setting "Apply Definition" to "False" using the palette and then reevaluate the input). As an example consider the 
product q[i,y]**p[i,z], it is not conserved and we get a messy answer: 
Commutator [q[i, y] **p[i, z], H] 

h 1 

i— **p[i, y] **p[i, z] -i (e 2 Z ft) ** **q[i, y] **q[i, z] 

m (q[i, x] 2 +q[i, y] 2 + q[i, z] 2 ) 372 

-» 

We next define the Runge Lenz vector A[i] with a formal definition: 



A[i] = L[i] xp[i] + z ® 2 m 



q[i] 




Using the NCcross function and symmetrizing in the two variables L and p, we work with 



A[i] = - NCcross [L[i] , p[i]] - - 



NCcross [ p [i] , L [i] ] + Z e 2 m 



q[i] 




We check a component of A to make sure it makes sense. 



A[i] [[3]] 



1 

- (-P[i# x] ** (-q[i, x] **p[i, z] ) +p[i, y] ** (-q[i, z] **p[i, y] ) - 
2 



p[i, x] **q[i, z] **p[i, x] +p[i, y] **q[i, y] **p[i, z] ) + 



1 



- (-(-q[i# x] **p[i, z]) **p[i, x] + (-q[i, z] **p[i, y] ) **p[i, y] + 
2 



e 2 m Z q [i , z ] 



q[i, y] **p[i, z] **p[i, y] -q[i, z] **p[i, x] **p[i, x] ) + 
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terml = Commutator [a [i] [ [3] ] , h] 

1 i 
2 



i le 2 Zi) * * p [i , z ] * * 



^/q[i, x] 2 +q[i, y] 2 +q[i, z] 2 



1 1 

— i (e 2 Zfi) ** **p[i, z] + 

y[q[i, x] 2 +q[i, y] 2 +q[i, z] 2 

1 

(e 2 Zfi 2 ) ** **q[i, z] - 

(q[i, x] 2 +q[i, y] 2 + q[i, z] 2 ) 372 

1 1 

— i (e 2 Zi) **p[i, x] ** **q[i, x] **q[i, z] 

2 " 

1 

— i (e 2 Zh) **p[i, y] ** **q[i, y] **q[i, z] 

2 



(q[if x] 


2 +qi 




y] 2 


+ q[i, z] 2 ) 3/2 








l 




(q[if x] 


2 +q 




y] 2 


+ q[i, z] 2 ) 372 








l 





i (e 2 Z h) **p[i, z] ** — - **q[i, z] **q[i, z] + 



1 

2 (q[if x]* +q[i, y]" +q[i, z 

1 1 



— i Z nj ** **p[i, xj **q[i, xj **q[i, zj + 

2 (q[i, x] 2 + q[i, y] 2 +q[i, z] 2 ) 372 

1 1 

— i (e 2 Zii) ** **p[i, y] **q[i, y] **q[i, z] - 

2 (q[i, x] 2 + q[i, y] 2 + q[i, z] 2 ) 372 



(q[i, x] 2 +q[i, y] 2 +q[i, z] 2 ) 372 



• P [if z] ** q [if x] ** q [i, x] 



(q[if x] 2 +q[i, y] 2 +q[i, z] 
1 1 



** p [if z] ** q [if y] ** q [if y] 



**p[i, z] **q[i, z] **q[i, z] 



2 (q[i, x] 2 + q[i, y] 2 +q[i, z] 2 ) 372 

With this expression, SimplifyQ does not produce a vanishing result since the expression on right locates p[i,z] 
differently in various terms and one needs to push this term to one side in ALL terms. Therefore at this stage it is better 
to push all the momentum type variables to the right extreme ( the left extreme works equally well), so that all 
operators are ordered in a similar way. After that we can set all non commutative products to simple products and 
simplify. The function pp and ppp achieve this ordering, and sim does the job of simplifying the expressions after 
replacing the **'s with simple *'s. 
pp[a_, m_] := PushOperatorRight [a, p[i, m] ] 

PPP[a_] •■= PPtPPtPPtaf x] , y] , z] 

sim[a_] := Simplify [ppp [a] /. {NonCommutativeMultiply -» Times} ] 
sim [terml] 



This can be examined more closely by examining the intermediate terms: 
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ppp[terml] /. {NonCommutativeMultiply -» Times} 



3e 2 Zft 2 q[i, x] 2 q[i, z] 3e 2 ZA 2 q[i,y] 2 q[i, z] 



2 (q[i, 


T 

x] 2 + q[i, y] 2 +q[i, z] 2 ) 572 


2 (q[i, x] 2 +q[i, y] 2 + q[i, z] 2 ) 572 




3 e 2 Zh 2 q[i, z] 3 


ie 2 Zip[i, z] q [i, x] 2 


2 (q[i 


, x] 2 + q[i, y] 2 + q[i, z] 2 ) 572 


(q[i, x] 2 +q[i, y] 2 +q[i, z] 2 ) 372 


i 


e 2 Z Jjp[i, z] q[i, y] 2 


3 e 2 Z h 2 q[i, z] 



(q[i, x] 2 + q[i, y] 2 + q[i, z] 2 ) 372 2 (q[i, x] 2 + q[i, y] 2 +q[i, z] 2 ) 372 



ie 2 Ziip[i,z]q[i,z] 2 ie 2 ZAp[i,z] 

+ 

(q[i, x] 2 +q[i, y] 2 +q[i, z] 2 ) 372 ^ q[if x] 2 +q[i# y]2 +q [i, z ] 2 

Simplify [%] 



The same scheme is demonstrated by pushing operators to the left (rather than right) here. 
ppl[a_, m_] := PushOperatorLef t [a, p[i, m] ] 
PPPl[a_] s= ppl[ppl[ppl[a, x] , y] , z] 

siml[a_] := Simplify [pppl [a] /. {NonCommutativeMultiply -» Times} ] 
siml [terml] 



■ (III) Fermions at work 

• For this Section Ensure that Fermionic Operators as well as Bra and Ket Vectors are Activated. 

Here we use 1 to represent spin up and - 1 to represent spin down. We now perform some test operations in the context 
of three problems using the Hubbard model Hamiltonian, given below. 
Clear[i, j , n, H] ; 

H = -,EZZ t[i ' j]*t[i, a]**f[j, a]+U^n f [i, l]**n f [i, -1]; 

i j a i 

. Problem 1 : [If, H] 

These examples require the user to input expressions carefully. For example, notice that we have a sum over <x in 
the first term of H above. We therefore cannot use an operator that has a spin specified by cr, because the package does 
not differentiate that one sigma is inside the sum and the other is outside. It will attempt to sum over both cr's. This is 
not so much a restriction as it is a warning against sloppy notation that could easily be understood by a human but will 
be misinterpreted by the package. The same restriction applies to site indices, such as i and j in this example. 
The example below demonstrates incorrect usage of indices; we are taking the commutator of a fermi operator with a 
site index that is identical to the indice of summation in the Hamiltonian. 
Commutator [ f [ i , ct] , H] 

^(05[1, a]) **n f [i, -1] **f[i, 1] - 

i 

^(U6[-l, a]) **ft[i, 1] **f[i, -1] **f[i, 1] ~X!ZZ t[i ' j] ** f ^' °] 

i i j a 

The above answer is not correct. 

Trying same the problem again using proper notation where the external indices are different from the ones used in H, 
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Commutator [f[k, al] , H] 

(U6[l, al]) **n f [k, -1] **f[k, 1] - 

(U6[-l, al]) **ft[k, 1] **f[k, -1] **f[k, 1] -^t[k, j] **f[j, al] 

j 

we get the correct answer. We can simplify the answer further by specifying a value for crl . 
Commutator [f[k, al] , H] / . crl -» 1 

U**n f [k, -1] **f[k, 1] -^Tt[k, j] **f[j, 1] 

j 

The best way of performing this calculation is just to start with the spin specified. In general this method is preferred 
because the package is then able to organize the operators using the value of the spins. It may therefore be necessary to 
reorder the output again if you specify the value of a spin after calculation. 
Commutator [f[k, 1], H] 

U**n f [k, -1] **f[k, 1] -^Tt[k, j] **f[j, 1] 

j 

Commutator [f[k, -1], H] 

-U** ft[k, 1] ** f [k, -1] ** f [k, 1] -^Tt[k, j] ** f [j, -1] 

j 

■ Problem 2 : [ff nf"", H] 

In this example we see that the answer is different in terms of operator ordering, depending on which method we use to 
perform the calculation. 

Commutator [f [k, al] ** n f [k, -al] , H] /. al -> 1 

U ** f t [k, -1] ** f [k, -1] ** f [k, 1] - 

U**ft[k, -1] **n f [k, 1] **f[k, -1] **f[k, 1] -^Tt[k, j] **n f [k, -1] **f[j, 1] + 

j 

^t[i, k] **ft[i, -1] **f[k, -1] **f[k, 1] -^Tt[k, j] **ft[k, -1] **f[j, -1] **f[k, 1] 

i j 

Commutator [f[k, 1] **n f [k, -1], H] 

U**n f [k, -1] **f[k, 1] -^Tt[k, j] **n f [k, -1] **f[j, 1] + 

j 

^t[i, k] **ft[i, -1] **f[k, -1] **f[k, 1] -^Tt[k, j] **ft[k, -1] **f[j, -1] **f[k, 1] 

i j 

The two results are not identical. However, if we reorder the operators in the top output we replicate the bottom output, 
which has the operators correctly ordered. 

StandardOrderQ [Commutator [f[k, al] **n f [k, -al], H] /. al -» 1] 

U**n f [k, -1] **f[k, 1] -^Tt[k, j] **n f [k, -1] **f[j, 1] + 

j 

^Tt[i, k] **ft[i, -1] **f[k, -1] **f[k, 1] -^Tt[k, j] **ft[k, -1] **f[j, -1] **f[k, 1] 

i j 

Performing the same calculation with a spin down particle we obtain the same answer with the spins reversed. 
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Commutator [f[k, -1] **n f [k, 1], H] 

-U**ft[k, 1] **f[k, -1] **f[k, 1] -^Tt[i, k] **ft[i, 1] **f[k, -1] **f[k, 1] + 

i 

^]t[k, j] **ft[k, 1] **f[j, -1] **f[k, 1] + ^Tt[k, j] **ft[k, 1] **f[k, -1] **f[j, 1] 

j j 



Problem 3: Hf + Jf + 2 | Vacuum} 



This demonstrates the properties of the Ket and Bra states, which need to be activated on the palette in addition to the 

Fermi, Bose or other relevant operators. 

We can perform the product directly. 

H8£t[l, 1] ** ft [2, -1] ** Ket [Vacuum] 

- (u | Vacuum) ) ** f t [2, -1] ** f t [1, 1] **n f [2, 1] - 
(u |vacuum)) **ft[l, -1] **ft[2, -1] **ft[l, 1] **f[l, -1] + 

^ ( | Vacuum) t[i, 1]) ** f t [2, -1] ** f t [i, 1] + ( [vacuum) t[i, 2]) ** f t [i, -1] ** f t [1, 1] - 

i i 

(U |vacuum)) **ft[2, -1] **ft[i, -1] **ft[l, 1] **ft[i, 1] **f[i, -1] **f[i, 1] + 

i 

XXX ( | Vacuum ) M 1 ' ^) ** f t [2, -1] ** f t [1, 1] ** f t [i, a] ** f [ j, a] 

i j a 

We see that the destruction operators annihilate the vacuum state 

f [i, cj] ** Ket [Vacuum] 



We could use the ProductQ operator here, or equivalently use the symbol ® between the operator and the ket. 

ProductQ[H, Ket [Vacuum]] 



Ket[l] = ft[l, 1] ** ft [2, -1] ** Ket [Vacuum] 

ft[l, 1] **ft[2, -1] ** [vacuum) 

ProductQ[H, ft[l, 1] ** ft [2, -1] ** Ket[Vacuum]] 

t [i, 1] ** f t [2, -1] ** f t [i, 1] ** | Vacuum) + t [i, 2] ** f t [i, -1] ** f t [1, 1] ** | Vacuum) 

i i 

ProductQ [H, Ket[l] ] 

[i, 1] ** f t [2, -1] ** f t [i, 1] ** | Vacuum) + [i, 2] ** f t [i, -1] ** f t [1, 1] ** | Vacuum) 

i i 

H®ft[l, 1] ** ft[2, -1] ** Ket[Vacuum] 

^t[i, 1] ** f t [2, -1] ** f t [i, 1] ** | Vacuum) + t [i, 2] ** f t [i, -1] ** f t [1, 1] ** | Vacuum) 

i i 

We thereby see that all destruction operators have been killed and the correct answer is obtained. 
Trying the above problem w/o ®, 
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H ** ft[l, 1] ** ft [2, -1] ** Ket[Vacuum] 

* * f t [ 1 , 1] * * f t [ 2 , -1] * * I Vacuum ) + 



uVft[i, 1] **f[i, 1] **ft[i, -1] **f[i, -1] 

i 

XZZ ft [i ' ct] ** f [j ' ct] t[i ' j] 



\ 1 

I 

**ft[l, 1] **ft[2, -1] ** I Vacuum) 



Note that we would not get the correct answer if we just used NonCommutativeMultiply instead 
of ® . Where it is appropriate to use NCM (**) rather than ((££)) might be confusing initially, 
but generally it is easy to deduce which one is appropriate. A rule of thumb is to useNCM in between single terms, 
and to use® in between larger expression involving many terms, Plus, Minus, or Sum functions). 

(IV A) Hubbard Model k - space Moment Calculations 
Clear[H]; H[l_, o_, p_, q_, r_] : = 

f+[1 ' CT] ** f t 1 ' CT 1 + (U/ Ns ) XXX fnP + (1 ' 1] ** f [P' ^ **ft[r-q, -1] **f[r, -1] 

la p q r 

This is the Hubbard Hamiltonian in k space, with Ns as the number of sites, and the labels p,q,l, r as the wavevectors. 
and £=e(k)- //, with /i as the chemical potential and band dispersion e(k). The dummy variables summed over in the 
right hand side (p,qj,l,cr) are declared on the left hand side explicitly. 

Throughout this notebook it is advisable to explicitly define the spin of the operators we are commuting with the 
Hamiltonian. This allows delta functions to be cancelled and thereby simplifies the result. Below we have computed 
the first moment using both spin up and spin down operators (1 for spin up, -1 for spin down). Also remember we 
must use new variables every time we commute with another Hamiltonian. 

We now compute the symmetric moments of the Hubbard model- these are defined 

as u> m = (-1)™ < {[H, [H, ... [H, /], ffl > with m nested commutators. Thus explicitly the first moment is the expecta- 
tion: oj\(k) = (-) < {[//, f(k, cr)], ft(k, <x)} >. Below we work out the operators that emerge from this nested commuta- 
tor and the final anticommutator. (Clearly the expectation of these requires the solution of the eigensystem that we do 
not address). 

A good reference for these moments is W. Nolting, Z. Physik 255, 25— 39 (1972), his Eqs(8-ll) correspond to our 
moments below after shifting oj m = M {m+l \ 

First Moment 

-Anticommutator [Commutator [H [1, a, c, q, r] , f[k, 1]], f t [k, 1]] 

> — ** n f [r, -1] + £[k] 
V NS 

Second Moment 
Anticommutator [ 

Commutator [H[l lr oi, x, &i, Pi], Commutator [H [1 , a, c, q, r] , f[k, 1]]], f t [k, 1]] 



> **n f [r, -1] + ) **n f [P!, -1] + ) ) **n £ [p lf -1] 



** f t [r + 9i, -1] ** f t [-9i + pi, -1] ** f [r, -1] ** f [p lf -1] 

r 6i Pl Ns 



ir 

V V V ** i t [-ei + p lt - l] ** f t [k + Si, l] ** f [r, - l] ** f [k - r + p lf 1] 

^^^NS 2 



r 0! Pi 

t2 



** f t [-q - 9i + Pi, -1] **ft[k + 9i, 1] **f[pi, -1] **f[k-q, 1] +f[k] 

q 9! Pi NS 
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Clearly a relabeling of variables would help combine terms here; it is left as an exercise for the reader. 
Third Moment 

The third moment can be calculated, but it requires further simplification. Due to the length of the result it is not shown 
below. Remove the semicolon at the end to see the large output. This calculation is long, and it takes around ten 
seconds to complete. 

- AntiCommutator [Commutator [H [A, J5, g, w, (.] , 

Commutator [H[l lf o lt x, 6i, Pi] , Commutator [H [1 , a, p, q, r] , f [k, 1] ] ] ] , ft [k, 1] ] ; 

(IV B) Hubbard Model r - space Moment Calculations 

Here the same calculation as above can be done in real space-finally taking Fourier transforms over the two fixed site 
indices rlj2 with wave vector k gives the moments 
co m (k). The results are equivalent to those in k space shown 

above. It is interesting to see these since the Pauli principle identities such as n ** n = 
n are utilized here and we can finally take Fourier transforms and sum over the fixed external sites (rl , r2 ) . 
Clear [H] ; 

H[j_] := (-1) * Z Z tk [J]'">[i] ft t k [j], u [jH **f[m[j], u[j]] - 

k[j] m[j] u[j] 

H ^ ^ ft[s[j], v[j]] **f[s[j], v[j]] - U ^ n f [t[j], 1] **n £ [t[j], -1] 

s[j] v[j] t[j] 

First Moment 

-AntiCommutator [Commutator [H [1] , f[rl, 1]], ft[r2, 1]] 

- (U 6 [rl, r2] ) ** n f [rl, - 1] - t r i, r2 - yii 6 [rl , r2] 

Second Moment 

AntiCommutator [Fold [Commutator, f[rl, 1], {H[l], H[2]}], ft[r2, 1]] 

(Ut r i, r2 ) **n f [rl, -1] + (Ut r i, r2 ) **n f [r2, -1] + (u 2 6[rl, r2 ] ) **n f [rl, -1] + 
2 (U/j6[rl, r2] ) **n f [rl, -1] + 2(it rl , r2 + ]T (Ut rl , m[2] 6[rl, r2] ) ** ft [rl, -1] ** f [m[2] , -1] - 

m[2] 

Z (Dtk[2],n <5[rl, r2] ) ** ft[k[2] , -1] ** f [rl, -1] + ]T t rl , m[1] t m[1] , r2 + u 2 <5[rl, r2] 

k[2] m[l] 

Third Moment 

(* Warning: Produces a large output *) 

-AntiCommutator [Fold [Commutator, f[rl, 1], {H [1] , H [2 ] , H [3] } ] , ft[r2, 1]]; 

(V) Spins at work: Eight Vertex Model Yang - Baxter Equation 
• Activate the Pauli matrix symbol 'V 

We use DiracQ to simplify some algebra that arises in the well known work of Baxter on the two dimensional 8 
vertex model (R. J. Baxter, Ann. Phys. vol 76, pi (1973)) . The transfer matrix of the model with a parameter u is 
written as: 

T (u) = Tr[L N>g (u) L N _ 1>g (u) ...L 1>g (u)] , (1) 

s 

where each of the objects L n g (u) is a scattering matrix 

a+b a-b c+d c-d 

Ln, g (u) = + + <o-g + oicfg. (2) 
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Here the Pauli matrices are defined for sites n=l,N and an extra "ghost" site "g". The Boltzmann weights a, b, c, and 
d are functions of the spectral parameter u. The fundamental commutation relation of the transfer matrices 
[T(u),T(v)]=0 defines the integrability of this model and also is the key to its solution. Baxter shows that this commut- 
ing transfer matrix relation is satisfied if the local relation 

L n>g2 [u] <g>L n>gl [u']®R glig2 [u"]= R gl , g2 [u n ] ®L n>g| [u'] ®L n>g2 [u] (Baxter's equation) 
is satisfied, where explicitly 

a + b a-b c + d c-d 

L„, g2 (u) = — + — o* b crl + —o* i2 <Tl + — oi 

a'+b 1 a'-b' c' + d 1 c'-d' 

L„, gl (U ') = — J— + o\ x (Tl + — J— (T X gi (Tl + — J— ( 3 ) 

a"+b" a"-b" c" + d" c"-d" 



+ 



oi. oi, + oi. crl + oi, ol . 



Let us use the functions of the DiracQ package to see if we can find the requirements on the u's such that (10) is 
satisfied. For convenience let us set the three distinct sites in the problem explicitly as: n=3 , g 2 = 1 , gi = 2 . This 
allows the code to automatically solve the Kronecker 8 functions, as shown below. For more information on use of the 
Kronecker 6 see the Tutorial. 
S[gi, g 2 ] 

<5[gi, 92] 

S[2, 1] 



n = 3;g 2 = l;g 1 = 2; 

a+ba-b c+d c-d 

Vg 2 ["] = + °l92, z] **a[n, z] + a[g 2 , x] **a[n, x] + a[g 2 , y] **a[n, y] ; 

2 2 2 2 

L n, gi [W] = 

a'+b'a'-b' c'+d' c'-d' 

+ a[gi, z] **a[n, z] + o[g lr x] **a[n, x] + °[gi, y] **a[n, y] ; 

2 2 2 2 

a''+b'' a''-b'' 

R gi,g 2 [W ' ' ] = + o[g lr z] ** a[g 2 , z] + 

2 2 

c''+d'' c''-d'' 

a[gi, x] **a[g 2 , x] + o[g lr y] **a[g 2 , y] ; 

2 2 

Let us first compute the product of the LHS Baxter's equation. For this we will use the function, which will apply the 
algebra of Pauli matrices and collect c terms. Because there are 8 terms in each scattering matrix the number of terms 
on both the right and left hand sides of equation is large (128 terms in all). Because of this we will suppress the output. 

LHS = L n , 92 [W] ® (L n , gi [W ] ®R gi ,g 2 [W " ] ) ; 

[W] 9 ["]); 
{Length [LHS] , Length [RHS] } 

{128, 128} 

We can take a peek at the output: 
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Short [LHS, 3] 

1 1 

— (a a' a") **a[l, z] **a[2, z] + — (a a' a") **a[l, z] **a[3, z] + 



4 4 

1 111 

- (a a' a") ** o[2 , z] ** a [3, z] + «170» + — d d' c" + — d c' d" + — c d' d" 
4 4 4 4 



Short [RHS, 3] 

1 1 

— (a a' a") **a[l, z] **a[2, z] + — (a a' a") **a[l, z] **a[3, z] + 



4 4 

1 111 

— (a a' a" ) ** a [2 , z] ** a [3 , z ] + «170» + — d d' c" + — d c' d" + — c d' d" 
4 4 4 4 



We now need to collect terms which involve the same set of operators. We can then set the sum of the coefficients of 
all the terms with identical operators to zero. We will now use the function TakeQPart, to extract the string of opera- 
tors from each term. By using the Union function we can also remove duplicate expressions. We are thereby left with a 

list of all the different strings of operators that appear anywhere throughout the expression. 
OperatorTerms = Union [TakeQPart [LHS - RHS] ] 

{a[l, x] **a[2, y] **a[3, z], a[l, x] **a[2, z] **a[3, y] , a[l, y] **a[2, x] **a[3, z] , 
a[l, y] **a[2, z] ** a[3, x] , a[l, z] **a[2, x] ** a[3, y] , a[l, z] **a[2, y] **a[3, x] } 

Length [OperatorTerms] 

6 

We can now use the QCoefficient function to find the coefficients of a string of operators. The use of this function is 
demonstrated below. 

QCoefficient [LHS - RHS, OperatorTerms [ [1] ] ] 

11111111 

— i d c' a" i c d' a" + — i c c' b" i d d' b" + — i b a' c" i a b' c" + — i a a' d" ibb' d" 

22222222 

The function has found all of the terms with operators that match the pattern of OperatorTerms [[1]] and summed their 
coefficients. This term must go to zero for the Yang - Baxter relation to hold. We can work this function into a loop 
and find all such expressions at once, as shown below. 

Do [expression [n] = QCoefficient [LHS - RHS, OperatorTerms [ [n] ] ] ; 
Print [expression [n] ==0], {n. Length [OperatorTerms] } ] 

11111111 

— i d c' a" i c d' a" + - i c c' b" i d d' b" + — i b a' c" i a b' c" + — i a a' d" ibb' d" == 

22222222 

11111111 

— i d a' a" + — i c b' a" i c a' b" + — idb'b" i b c' c" + — iad'c" i a c' d" + — i b d' d" == 

22222222 

11111111 

— i d c' a" i c d' a" i c c' b" + — i d d' b" i b a' c" + — i a b' c" + — i a a' d" ibb' d" == 

22222222 

11111111 

— i d a' a" i c b' a" + — i c a' b" + — i d b' b" + — i b c' c" + — i a d' c" i a c' d" i b d' d" == 

22222222 

11111111 

— i a c' a" i b d' a" i b c' b" + — i a d' b" i c a' c" + — i d b' c" + — i d a' d" i c b' d" == 

22222222 

11111111 

— i a c' a" i b d' a" + — i b c' b" + — i a d' b" + — i c a' c" + — idb'c" i d a' d" i c b' d" == 

22222222 

These are the final equations. By recombining them and simplifying we reproduce Baxter's equations. 
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Expand [ (expression [2] + expression [4] ) / i] == 
Expand [ (expression [ 1] + expression [3 ] ) / i] ==0 
Expand [ (expression [5] + expression [6] ) / i] == 
Expand [ (expression [2] - expression [4] ) / i] == 
Expand [ (expression [ 1] - expression [3 ] ) / i] == 
Expand [ (expression [5] - expression [6] ) / i] == 

-d a' a" + d b' b" + a d' c" - a c' d" == 
d c' a" - c d' a" + a a' d" - b b' d" == 
-b d' a" + a d' b" + db'c"-c b' d" == 
c b' a" - c a' b" - b c' c" + b d' d" == 
c c' b" - d d' b" + b a' c" - a b' c" == 
a c' a" - b c' b" - c a' c" + d a' d" == 

These are six homogeneous equations relating the twelve parameters (a,b,c,d) etc. Baxter proceeds to give an elegant 
analysis of these by parameterization in terms of the spectral parameters u, u', u". The solution is well described in 
the original paper; we see above that DiracQ helps in automating the tedious calculations. 
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■ (VI) ID Hubbard Model SSS Matrix Problem 
• Activate the Pauli matrix symbol "<r" 

We next demonstrate the use of DiracQ in the context of the 1-d Hubbard model. It satisfies the Yang Baxter relation, 
analogous to the 8 vertex model discussed in the previous example, and the proof of this is given in B. S. Shastry, 
Phys. Rev. Letts, vol 56, 1529, 2453 (1986). A relatively simple analytical proof is also to be found in B. S. Shastry, 
J. Stat. Phys. vol 50, 57 (1988), using the decorated star triangle relations. We recall this result before proceeding 
further. It may be noted that the R matrix of this model, defined below, is currently of great interest in connection with 
the problem of the ADS/CFT correspondence for N = 4 SUSY Yang-Mills theory in four dimensions following the 
work of N. Beisert, J.Stat.Mech. 0701 (2007) P01017, arXiv:nlin/0610017. 

The 1 -dimensional Hubbard model Hamiltonian using the Jordan Wigner transformation and with suitable boundary 
conditions can be written in terms of Pauli spin matrices as 

^ 1 

H = 2j Hn+ l,n ; H n+1>n = (0- n + 0- n+1 +0- n+1 + 0- n ) + (T n + T n+1 +T n+1 + T n ) + " UCT^T^. (4) 
n 

This corresponds to the quantum Hamiltonian, commuting with the transfer matrix of two Free fermi models coupled 
in a specific way. We start with a single Pauli species scattering matrix : 

a + b a-b c , 

if (6) = — + — rfo-f + - {tn'trf + tftrf) (5) 

Here the parameters a b c satisfy the Free Fermi condition a 2 +b 2 = c 2 . This is conveniently parametrized by setting 

a = c Cos (6) b = c Sin (ff) . For the two sites 1 and 2 we define 

a^Cosffli], b^Sinf^], a 2 =Cos[0 2 ], b 2 = Sin[0 2 ] and denote (Wj = (0j, hj)). 

We may write two sets of scattering matrices ff) and ffj for the two species at the same 6 : 

Lj^lff (0)<g>l|?(0) (6) 
The scattering operator for the Hubbard model is constucted as 

L3„(W„) = 1 3 „(0„)I„ (h n ), 

The transfer matrix can now be written as 

T (W) = Tr[L No (W) L N _, (W) ... L, (W)] . (7) 

o 

With these definitions, the transfer matrix forms a commuting family [T (Wj), T(W 2 )] = 
provided the scattering matrix is proven to satisfy the Yang - Baxter equation (W n = (0„, h n )) 

L31 (Wi) L 32 (W 2 ) Su (W 2 I Wi) = S12 (W 2 I WOL32 (W 2 ) L31 (Wi). (8) 

The S operator found by Shastry is given by: 

S12 (W 2 I Wi) = p[cos (6 2 + 6»i) cos (h 2 - hi) 1 12 (0 2 - Oi) + cos (0 2 - 6»i) sinh (h 2 - h x ) \ n (6 2 + Q\)o\ t|]. (9) 

An extra and crucial sufficiency condition is that the fields h„ (n = 1 , 2) must satisfy 

Sinh[2h n ] U 

= - , (10) 

a n b n 2 

where U is the interaction term in the Hubbard Hamiltonian. The peculiarity of 
this model is that S (W 2 | Wj ) depends on the spectral parameters Wj and W 2 separately, 

and cannot be written in terms of the difference of the spectral parameters. This relation 
was proved by Shastry analytically using the decorated star triangle relations. 
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A more (technically) ambitious question is to ask is if we 
can construct a more general inhomogeneous model than the Hubbard model, 
with a transfer matrix T (W | Wj , ... W„) depending on the N fixed parameters Wj 

T = Tr[S No (W, | W 2 ) S N _, (W | W N _,) ... S, (W | WO]. 



(11) 







so that on setting Wj to 0, we get back the Hubbard transfer matrix 
above. The inhomogeneous matrix T commutes for different weights according to 

[T (W | Wj ... W N ), T(W' | (Wj ... W N )] = 0. 

The sufficiency condition for this is the more general Yang Baxter relation 

S31 (W, | W 3 ) S 32 (W 2 | W 3 ) Su (W 2 | Wi) = Su (W 2 | Wi) S 32 (W 2 | W 3 ) S 31 (W, I W 3 ). 



(12) 



(13) 



This is Eq( 4.14) of the paper and several typical matrix elements were verified by Shastry by explicit calculation, who 
conjectured that this equation is true for all matrix elements. Since the number of terms is very large (see below), an 
exhaustive enumerative proof was not given. We will provide such an enumerative proof using DiracQ next. 

For nonzero Wj, an analytical proof was actually found much later in the work of M. Shiroishi and M. Wadati - 
Journal of the Physical Society of Japan, Vol 64, 57-63 (1995), and thus the following enumerative proof is not the 
first proof. However, it does demonstrate the power of DiracQ quite well. The large number of terms involved make 
the algebra prohibitive by a hand calculation for most of us! Using the package however we are able to compute the 
product of three S operators in a number of minutes. 

The first step is to input the form for the scattering matrix and for the S operators. Here, we denote the two separate 
species of Pauli matrices not by use of another symbol, but rather by using a third index to denote the species of the 
operator, as shown below. 



In this notation the scattering matrix and the S operator are written below as general functions. 

Clear [A, 1, S] ; 
l u _,v_[n_, p_, q_] : = 

(A[n, 1, p, q] +A[n, -1, p, q] ct[v, z, 1] **a[u, z, 1] +a[v, x, 1] **a[u, x, 1] + 

<?[v, y, 1] ** ct[u, y, 1] ) ® (A[n, 1, p, q] + A[n, -1, p, q] ct[v, z, 2] ** a[u, z, 2] + 
a[v, x, 2] **ct[u, x, 2] +cr[v, y, 2] **a[u, y, 2]); 
S u _,v_ : = Cos [6 V + 6 U ] Cosh [h v - h u ] l u , v [- 1» 9 V , U ] + 

Cos[e v -6 u ] Sinh[h v -h u ] l UiV [l, e v , 6 U ] ct[v, z, 1] a[v, z, 2] 

The Boltzmann weights in the scattering matrix are not written explicitly because this increases the number of terms 
involved in the product and their definition will not be relevant until after the S operator products have been calcu- 
lated. They are therefore represented by the function A, which will be defined after the product of the three S operators 
on each side of (36) has been computed. We will set 
rulesimplify = { h ± _ -> 1 / 2 ArcSinh [U / 2 Cos [e ± ] Sin [e±] ] , , 

A[n_, m_, p_, q_] ->Cos[p + nq] + m Sin [p + n q] } 
after the calculation with the operators is performed, and it is easy to see that this replacement reproduces the scatter- 
ing operator above with a scale factor of 2 which cancels out. Here we compute the product of three S , matrices. Each 
S matrix contains about thirty terms, we we wexpect the product to contain on the order of 27,000 terms. This expres- 
sion is so large as to be prohibitive to compute by hand, 
{time, LHS } = Timing [S 3 ,i ® (S 3 , 2 ® Si, 2 ) ] ; 
time / 60 



^-»<r[l,z, 1], t?-»o-[1,z,2] 



(14) 



2.79112 
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The timing function informs us that this operation takes around three minutes on our testing machine. 
Length [LHS] 
27 360 

We see that there are approximately 32,000 terms in the final product. Below an arbitrary example of a single term is 
shown. We expect similar timing and length for the RHS. 

LHS [ [6] ] 

a[-i, i, ei, e 3 ] 2 A[-i, i, e 2 , e!] 2 A[-i, 1, e 2 , e 3 ] 2 

Cos [Oi + 2 ] Cos [9x + 3 ] Cos [e 2 + 3 ] Cosh [h! - h 2 ] Cosh [hj - h 3 ] Cosh [h 2 - h 3 ] 
We now compute the RHS 
RHS = Si, 2 ® (S 3 , 2 ® S 3 ,i) ; 
Length [RHS] 
27 360 

We have now computed both the left and right hand sides of equation 13. If the equation holds, they are equal. We 
verify this by subtracting the two and returnign zero, 
subtraction = LHS - RHS; 

Length [subtraction] 

13 472 

We would have hoped that when we subtracted the LHS from the RHS that the result was zero. It appears we will have 
to do more manipulation to show that the LHS is equal to the RHS. We now have a very large expression, and we need 
to manipulate it in several ways. We need to find all distinct strings of operators and then find the coefficient of all of 
each string. We must then show that every one of these coefficients vanishes. DiracQ contains functions that will 
perform the first two. However, were we to use the functions outright the computation time would be very long. Every 
expression input into a DiracQ function is first organized into lists that can be manipulated (see the ' Brief Explanation 
of Form' section). For a large expression such as the one we are working with, organizing the expression takes a very 
long time. The expression must then be recompiled into and understandable form. We can circumvent some of these 
steps by organizing an expression once and specifiying to subsequent functions that the input has already been orga- 
nized and that this step should be skipped. This is done by specifying' Organized Expression -> True' as an option. 
This process is demonstrated below. 
? Organize 



Organize is the function that enables the DiracQ package to understand user input. Organize takes a mathematical 
expression as input and yields a nested list that contains the atoms of the input ordered according to 
their properties. Numbers, summed indices, c numbers, and q numbers are separated into groups. Each 
term of the input separated by plus sign constitutes a separate list of items in the output. Example: 

Organize[(tt) V (c tt)*(q tt)J 

index 

={{B,{index|,{c tt),(q it}}) 

Organize^) ]T (c tt]>(q tt 1 )+(tt 2 ) £ (c tt 2 )*(q tt 2 )] 

indexi index 2 

Hltti.lindeXiUc tt^lq tt 1 )),{tt 2 ,{index 2 ),{c tt 2 ),{q tt 2 )U 

For a more in depth explanation see the DiracQ writeup notebook. 

{time, organizedsubtraction} = Timing [Organize [subtraction] ] ; 
time / 60 

1.37672 
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We see that on our test machine the Organize function takes a few minutes to evaluate and is computationally quite 
expensive. Below we use the TakeQPart function to find the strings of operators that appear in OrganizedSubtraction, 
specifying that the input is already an organized expression to save some computational time. 
Operator-Terms = Union [TakeQPart [organizedsubtraction, OrganizedExpression -* True] ] ; 

Length [OperatorTerms] 

472 

OperatorTerms [ [300] ] 

a[l, y, 1] **a[2, z, 1] **a[3, x, 1] **o[2, z, 2] **a[3, z, 2] 

The above output shows that there are 472 distinct strings of operators found in the expression and an example of one 
such string. We can now use the QCoefficient function to find the coefficients of a string of operators. 
? QCoefficient 



QCoeffficient[expression,form] will scan the expression for terms containing a string of operators that match 
'form'. The function will output the coefficient of the operator(s) specified by 'form'. Only exact 
matches of the string of operators are found. If several terms are found in expression containing 
terms that match 'form' the output will be a sum of the coefficients of the specified operators. 

The result is not printed due to the length and complexity, though we do show the first term in the expression. See the 
appendix section, or remove the semicolon at the end to see the expression. 
sum[l] = 

QCoefficient [organizedsubtraction, OperatorTerms [ [1] ] , OrganizedExpression -* True] ; 
Length [sum[l] ] 

24 

sum[l] [[1]] 

4 a [ - 1 , - 1 , e 2 , e 3 ] a [ i , - 1, ei, e 3 ] cos [0! - e 2 ] 

Cos [0i - 3 ] Cos [9 2 +63] Cosh [h 2 - h 3 ] Sinh [hj - h 2 ] Sinh [hi - h 3 ] 
We now apply the definition of the Boltzmann weights. 

rulesimplify = 

{ hi_ -> 1 / 2 ArcSinh [U / 2 Cos [6i] Sin [0i] ] , A[n_, m_, p_, q_] -> Cos [p + n q] + m Sin [p + n q] J ; 
Simplify [sum [1] /. rulesimplify] 



We see that the coefficients of this one string of operators simplify to zero. We must show that each such term van- 
ishes for all 472 strings of operators. The command below will evaluate each such term. 

Do[sum[i] = Simplify [QCoefficient [organizedsubtraction, 

OperatorTerms [ [i] ] , OrganizedExpression -» True] /. rulesimplify], {i, 1, 472}] 



Examplebook_VO.nb | 17 



Table[sum[i] , {i, 1, 472}] 
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o, 


o, 


o, 


o; 





Thus we see that every term simplifies to zero and thus proves the result. 



■ (VII) Hubbard Model Currents 

Below we define the Hamiltonian of the Hubbard model using fermionic operators. The system is taken to have 
periodic boundary conditions such that (N + a) a where N is the number of sites. We also write two currents 
which have been shown to commute with the Hamiltonian in the papers B.S.Shastry Phys. Rev. Letts, vol 56, 
1529, 2453 (1986) and especially J.Stat.Phys, vol 50, 

57 (1988). Here it is shown that for a system with odd number of sites we have a conservation law termed j A 
(defined below) that stands apart from the other conservation laws. A current j B is also found (defined below) 
and unlike j A generalizes to many more currents since it is the first non trivial member of the currents 
contained in the transfer matrix. Note that j A corresponds to infinite ranged hopping of the doubly occupied 
sites adding to short ranged hops of fermions, whereas j B corresponds to second neighbor hopping. 

We will now verify the commutation of j A for odd number of sites using DiracQ and the 
commutation of j'b for both odd and even sites. Notice that we write the Hamiltonian as well as the 
currents as functions of the number of sites NN (N is a predefined symbol in Mathematicd) . This 
will allow us to vary the number of sites between even and odd values. 

Clear [H] ; Remove [ j ] ; 

HH 

H[NN_] := -t ^Sum[ (f t [n + 1, a] ** f [n, a] + ft [n, a] ** f [n + 1, a] ) , {a, - 1, 1, 2}] + 

n=l 

HH 

U^n f [m, 1] **n £ [m, -1] 

m=l 

HH 

j A [NN_] := it ^]sum[ (-ft [1, a] ** f [1 + 1, a] + ft [1 + 1, a] ** f [1, a] ) , {a, -1, 1, 2}] + 

1=1 

HH HH-1 

iU^^ (-l) 1 ft[l + r, 1] **ft[l + r, -1] **f [r, -1] **f [r, 1]; 

r=l 1=1 

HH 

j B [NN_] ; = i t ^ Sum [ft [o + 2, a] **f[o, a] - ft[o, a] **f[o + 2, a], {a, -1, 1, 2}] + 

o=l 

HH 

ill ^ Sum [ft [o + 1, a] **f[o, a] -ft[o, a] **f[o + l, a], {a, -1, 1, 2}] + 

o=l 
HH 

iU ^]sum[ (f t [o, a] ** (f [o + 1, a] - f [o - 1, a] ) - (ft [o + 1, a] - ft [o - 1, a]) ** f [o, a] ) 

o=l 

nf [o, -o] , {a, -1, 1, 2}] ; 

First we will perform the commutator of the Hamiltonian with j A for a system 

with 5 sites. As the number of sites increases so does the number of operations performed, 
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and therefore we will only perform these commutators for relatively small systems. Even so the 
computing time is nonnegligible. We must include the periodic boundary conditions of our 
system. These are written as replacement rules. The replacement rules must be specified to each term, 

and cannot be specified at the end of computation, as this will lead to an erroneous answer. 

PBCrule [n_] : = {f [i_, a_] /; i > n -» f[i-n, a], f[i_, a_] /;i<l-»f[i + n, a], 
ft [i_, a_] /; i > n -» ft [i - n, a] , ft [i_, a_] /; i < 1 -» ft [i + n, a] } 

Timing [Commutator [H [5] /. PBCrule [5], j a [5] /. PBCrule[5]]] 

{20.779, 0} 

Performing the same computation with a four site sytem we see that the result is nonvanishing, as we expected. 
Commutator [H [4] /. PBCrule[4], j A [4] /. PBCrule[4]] 

2 i (tU) ** f t [1, -1] ** f t [1, 1] ** f [1, -1] ** f [2, 1] + 
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(tU) 


** ft [1, 


-1] 


** f t [1, 


1] 


* * f [ 2 , 


-1] 


** f [1, 


1] 
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(to) 


** ft [1, 


-1] 


** f t [2, 


1] 


* * f [ 2 , 


-1] 


* * f [ 2 , 


1] 
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(t U) 


** ft [1, 


-1] 


** f t [4 , 


i: 


** f [1, 


-1] 


** f [1, 


i: 




2 
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(t U) 


** f t [2, 


-1] 


** f t [1, 


i; 


* * f [ 2 , 
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** f [2, 
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(t U) 


** f t [2, 
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** f t [2, 


i] 


** f [2, 


-1] 


* * f [ 3 , 


i] 
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(t U) 


** f t [2, 
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** f t [2, 
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* * f [ 3 , 
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** f [2, 


i] 
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(t U) 


** f t [2, 
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** ft [3, 
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* * f [ 3 , 
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(tU) 


** f t [3, 
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** ft [2, 
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* * f [3 , 


-1] 
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i] 
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(to) 


** f t [3, 


-1] 


** ft [3, 


ij 


** f [3, 


-1] 


** f [4, 


i] 
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(to) 


** f t [3, 


-1] 


** ft [3, 


i] 


** f [4, 


-1] 


** f [ 3 , 


i] 
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(to) 


** f t [3, 
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** f t [4 , 
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i] 
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** f t [4, 
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** f t [1, 
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** f [1, 


i] 
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** f t [3, 


i] 


** f [4, 


-1] 
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(t U) 


** f t [4, 
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** ft [4, 


i] 


** f [1, 


-1] 


** f [4, 


i] 
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2 
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(to) 


** f t [4, 


-1] 


** f t [4 , 


i] 


** f [4, 


-1] 


** f [1, 


i] 





Again, a system with an odd number of sites commutes with the Hamiltonian. 
Timing [Commutator [H [7] /. PBCrule [7], j a [7] /. PBCrule[7]]] 

{59.062, 0} 

We now verify that j B commutes with the Hamiltonian for both odd and even number of sites . 
Timing [Commutator [H [5] /. PBCrule [5], j B [5] /. PBCrule [5]]] 

{42.806, 0} 

Timing [Commutator [H [6] /. PBCrule [6], j B [6] /. PBCrule [6]]] 

{62.323, 0} 

Here we investigate whether the two currents commute with each other. The results suggest that for odd N the two 

currents commute but that they do not commute for even N. 

Timing [Commutator [j A [5] /. PBCrule[5], jB[5] /. PBCrule[5]]] 

{95.364, 0} 

From the computation below we see that the commutator is non zero for even sites, 
evensitecommutator = Commutator [j a [4] /. PBCrule[4], j B [4] /. PBCrule[4]]; 
Short [evens itecommut at or , 2] 

2 (tU) **ft[l, -1] **ft[l, 1] **f[l, -1] **f[3, 1] - 
2 «l» + «44»+2 (tU) ** f t [4, -1] ** f t [4, 1] ** f [4, -1] ** f [2, 1] 

In view of the above results, it seems very likely that [j A , T] = for odd number of sites, where T is the transfer matrix 
of the Hubbard model in the Fermi representation, whereas the commutator is non zero for even sites. 

■ (VIII) Construction of Cluster Hamiltonians of the Hubbard model using Bras and Kets 
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We illustrate the use of the Bra[ Vacuum] and Ket[ Vacuum] symbols, for creating the space of allowed states in the 
Hubbard model on a small cluster, and furhter to operate the Hamiltonian on these to generate its matrix representa- 
tion. Let us consider a 4 site cluster defined in the diagram, with the Hubbard Hamiltonian hopping parameters (tl, t2 , 
U). For illustration we confine ourselves to the sector with Sz=0 and 4 particles, i.e. the half filled state, where the 
number of basis states is 36. 



(1) tl (2) 

tl t2 tl 

(3) tl (4) 

4-Site Hubbard Cluster with nearest neighbour hops tl , and second neighbour hops t2. 

in[3]:= hop[i_, j_] := - Sum[f t [i, a] ** f [j , a] + ft [j , a] ** f [i, a] , {a, - 1, 1, 2}] 
in[4]:= pot[i_] := Un f [i, 1] **n f [i, -1] 

in[5]:= Hcluster = tl (hop[l, 2] +hop[2, 3] +hop[3, 4] +hop[4, 1]) + 

t2 (hop[l, 4] +hop[2, 3]) + pot[l] +pot[2] +pot[3] +pot[4] ; 

We now define the creation and destruction operators for pairs of sites with a given spin projection, 

in[6]:= twofermi[i_, j_, a_] = ft[i, a] ** f t [ j , a]; 
twof ermidest [i_, j_, a_] = f[j, a] **f[i, a]; 

Let us see the explicit form of a basis ket. The function StandardOrderQ organizes a Fermi operator product, its 
convention is to write these with increasing site labels from left to right, and puts the down spin blocks to the left of 
the up spin blocks. 

in[8]:= StandardOrderQ [twof ermi [ 1, 2, 1] ** twof ermi [2 , 3, -1] ** Ket [Vacuum] ] 

Out[8]= ft[2, -1] **ft[3, -1] **ft[l, 1] **ft[2, 1] ** |vacuum) 

We next produce the 36 basis kets {¥[1], .., ¥[36]} and their adjoint 36 basis bras {¥B[1], .., <PB[36]} 

in[9]:= ii = 1; Do[5[ii] = StandardOrderQ [twof ermi [i, j, 1] ** twof ermi [k, 1, -1] ** Ket [Vacuum] ] ; 
ii = ii + 1, {i, 1, 4}, {j, 1, i-1}, {k, 1, 4}, {1, 1, k-1}]; 
ii = 1; Do[5B[ii] = StandardOrderQ [ 

Bra [Vacuum] ** twof ermidest [i , j, 1] ** twof ermidest [k, 1, -1]]; 
ii = ii + 1, {i, 1, 4}, {j, 1, i-1}, {k, 1, 4}, {1, 1, k-1}]; 

It is useful to inspect a typical basis state pair 

ln[11]:= SB [9] 

Out[ii]= ft[2, -1] **ft[3, -1] **ft[l, 1] **ft[3, 1] ** [vacuum) 

ln[12]:= 9iB [18] 

Out[i2]= (vacuum | **f[3, -1] **f [4, -1] **f[2, 1] **f[3, 1] 

We next illustrate the action of the Hubbard Hamiltonian on these states, here the function StandardOrderQ does the 
job of simplifying expressions bu pushing the destruction operators to the extreme right and hitting the Ket[Vacuum], 
which is annihilated, or a similar story for the Bra[Vacuum] that is annihilated by creation operators. A brief glance at 
the output shows that we get back the basis states multiplied by coefficients that depend on tl,t2 and U. 
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in[i3]:= Short [StandardOrderQ[Hcluster ** S[9] ] , 3] 

Out[i3]//Short= -tl ** f t [ 1, -1] ** f t [3 , -1] ** f t [ 1, 1] ** f t [3 , 1] ** | Vacuum) - 

tl ** f t [2, -1] ** f t [3, -1] ** f t [1, 1] ** f t [2, 1] ** | Vacuum) - «1» - «1» + «1» - 
«1» - «1» + «1» + U ** f t [2, -1] ** f t [3, -1] ** f t [1, 1] ** f t [3, 1] ** | Vacuum) 
in[i4]:= Short [StandardOrderQ[$B[ 18] ** Hcluster] , 2] 

Out[i4]//Short= 1 1 * * ( Vacuum | * * f [ 1 , -1] * * f [ 3 , -1] * * f [ 2 , 1] * * f [ 3 , 1] - « 1 » - « 1 » - 

« 1 » + « 1 » - « 1 » + U * * ( Vacuum | * * f [ 3 , -1] * * f [ 4 , -1] * * f [ 2 , 1] * * f [ 3 , 1] 

At this point, there are two ways to construct the Hamiltonian matrix. Starting with expressions of the type StandardOr- 
derQ[Hcluster***P[i]], we can take the inner product with the bras ¥B[j]. As an alternative, we can pick off the 
coefficients using the QCoefficient function, this is faster and hence prefereable in most situatiions. We illustrate both 
methods and display the time saved by the second method. First we need a rule that collapses the fundamental inner- 
product to unity: 

in[i5]:= ruleinnerproduct = {Bra [Vacuum] ** Ket [Vacuum] -» 1}; 

in[i6]:= Timing [Do [res [i] = StandardOrderQ [Hcluster ** & [i] ] , {i, 1, 36}]] 

Out[i6]= {4.5608, Null} 

The array res[i] stores the resultant of the action of the Hamiltonian on the ith ket. 

The two methods to creat the Hamiltonian matrix use the innerproduct with the *PB[j], i.e. the Bra vector, or picking 
out the coefficients of *P[j] in res[i]. These are called matrixl and matrix2 respectively. 

matrixl[j_, i_] := StandardOrderQ [SB [j ] **res[i]] /. ruleinnerproduct 
matrix2[j_, i_] := QCoefficient [res [i] , 5[j]] 

ln[34]:= 

Rl = Timing [Table [matrixl [k, i], {k, 1, 36}, {i, 1, 36}]]; 
R2 = Timing [Table [matrix2[k, i] , {k, 1, 36}, {i, 1, 36}]]; 

ln[36]:= Rl[[l]] 

Out[36]= 401.36 

ln[37]:= R2[[l]] 

Out[37]= 7.00216 

in[39]:= Short [Rl [ [2 ] ] -R2[[2]], 3] 

Out[39]//Short= {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, «34», {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}} 

We see that the two methods give the same resulting matrix, but the second method is considerably faster. We next 
define the ClusterHamiltonian and check its eigenvalues in some simple cases. For larger systems, it is clearly more 
advantageous to use the sparse matrix notation in such problems, but we will not pursue it here. 
in[40]:= ClusterHamiltonian [tl_, t2_, U_] = Rl[[2]]; 

in[4i]:= Eigenvalues [ClusterHamiltonian [1, 0, 0]] 

Out[4i]= {-4, -4, -4, -4, 4, 4, 4, 4, -2, -2, -2, -2, -2, -2, 

-2, -2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} 

in[42]:= Eigenvalues [ClusterHamiltonian [ , 0, U] ] 

Out[42]= {0, 0, 0, 0, 0, 0, U, U, U, U, U, U, U, U, U, U, U, 

U, U, U, U, U, U, U, U, U, U, U, U, U, 2U, 2U, 2 U, 2 U, 2 U, 2 U} 

The first set of eigenvalues correspond to a non interacting model with only nearest neighbour hoppings, and is easily 
seen to be correct. The second set of eigenvalues correspond to the local limit, where we have 6 states in the U=0 
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sector (Lowest Hubbard Band), 6 states in the U=2 sector (Second Hubbard band) and the rest with U=l, i.e. the First 
Hubbard Band. 



Appendices 

■ (I) Explanation of Form 

The strength of the system lies in the ability to organize, rearrange, and manipulate terms in an expression based on 
their properties. The method of organization provides a simple language for the manipulation of terms. An understand- 
ing of the language used by the code is not essential for most operations. The user is only required to input algebraic 
expressions in a logical and minimally restricted form to effectively use the functions of this package. The language of 
organization will briefly be covered here for posterity. 

Mathematica organizes all input in terms of lists. Each list has a head. The head of a list is the function that is to be 
applied to the list. The example below shows how Mathematica organizes some algebraic input into lists and functions. 
FullForm[5 £c[i] b[i] ** h\ [j] + 2 £d[i] bf [i] ** b[j]] 

i i 

Plus[Times[5, Sum[Times[cLi], NonCommutativeMultiply[bLi], b\[Dagger][j]]], i]], 
Times[2, Sum[Times|dLi], NonCommutativeMultiply[b\[Dagger]Li], b[j]]], i]]] 

It is the job of the Organize function to scan through this list and extract the relevant information. The Organize 
function is therefore at the heart of every operation of the code. This function reads input such as the example above 
and organizes it based on the properties of the terms encountered. Information encountered by the Organize function is 
stored in nested lists. Different types of objects are placed in different places within the larger list created by the 
Organize function. Organize recognizes four types of objects : numbers, sums, constants, and operators. Each of these 
types of objects is placed in a separate nested list. The easiest way to understand this ordering is through example 
Organize^ ^c[i] b[i] ** bt [j]] 

i 

{{5,{i),{c[i]),{b[i], btUJlll 



{{5, {i}, {c[i]>, {b[i], b|0]}}} 



T 

Numbers 



r t 

Constants Operators 



Sunnec 
Indices 



Fig. (1): The output of organize is a list of lists. Each list corresponds to a specific type of term. 



Organize^ ^]c[i] b[i] **bt[j] +2^]d[i] bt[i] **b[j]] 

i i 

{{5, {i}, {c[i]}, {b[i] , bt[j] }}, {2, {i}, {d[i]}, {bt[i] , b[j] }}} 



U5, {i}, {c[i]}, {b[i], btO]}}, {2, 0, {d[i]}, {bt[i], bO]}}} 

T T 

First Term Second Term 
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Fig. (2) : When the input of Organize is the addition of several terms, 
each term is contained in a list within the larger list. 

When a function is encountered that cannot be decomposed further, the Organize function determines whether this 
function contains any operators. If so, the package uses a special notation to signify that it has found a function of 
operators that cannot be decomposed further. The notation function[a, {b}] is used, where a is the function that cannot 
be decomposed, and b is a list of the operators on which the function depends. As it currently stands not all functions 
can be read and placed in this notation, but any function involving operators to different powers or exponential 
functions of operators can be understood. Using this format the user is free to define commutators of more complicated 
functions of operators as required. The notation used is identical to that used to add an operator, with the caveat that 
the definition will be written using the "function" notation described above. An example of an organized function such 
as described above is given below. 
Organize [a e t[i] iW+qlj]] 

{{1, {}, {a}, { function [ e q[ jl , {q[j]}], function [e q[il t[i] , {q[i]}]}}} 

Humanize is the functional opposite of organize in that Humanize reads the language created by the Organize function and 
recreates a mathematical form that a user can understand. This function is demonstrated below. 
Humanize[{{5, {i}, {c[i]}, {b[i], bt[j]}}, {2, {i}, {d[i]}, {bt [i] , b[j]}}}] 

5^Tc[i] **b[i] **bt[j] +2^d[i] **bt[i] **b[j] 

i i 

With this simple language the manipulation of certain types of terms can be performed with greater ease. Rather 
than manipulating operators as we come across them, we collect every operator in the expression. We can then per- 
form algorithms on the lists of operators to combine and manipulate them as necessary. It is necessary to collect 
numbers and constants separately so that numbers can be combined as necessary and constants, which may depend on 
indices of summation, can be placed inside a sum. Also, the ability to identify summed indices allows us to evaluate 
delta functions. Every function in the package uses this form to organize input. 
■ (II) Demonstration Problem Output (Produces large outputs hence not meant for printing) 



